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Statistical DNA models available in the literature are often effective models where the base-pair 
state only (unbroken or broken) is considered. Because of a decrease by a factor of 30 of the effective 
bending rigidity of a sequence of broken bonds, or bubble, compared to the double stranded state, the 
inclusion of the molecular conformational degrees of freedom in a more general mesoscopic model is 
needed. In this paper we do so by presenting a ID Ising model, which describes the internal base pair 
states, coupled to a discrete worm like chain model describing the chain configurations [J. Palmeri, 
M. Manghi, and N. Destainville, Phys. Rev. Lett. 99, 088103 (2007)]. This coupled model is 
exactly solved using a transfer matrix technique that presents an analogy with the path integral 
treatment of a quantum two-state diatomic molecule. When the chain fluctuations are integrated 
out, the denaturation transition temperature and width emerge naturally as an explicit function 
of the model parameters of a well defined Hamiltonian, revealing that the transition is driven by 
the difference in bending (entropy dominated) free energy between bubble and double-stranded 
segments. The calculated melting curve (fraction of open base pairs) is in good agreement with the 
experimental melting profile of polydA-polydT and, by inserting the experimentally known bending 
rigidities, leads to physically reasonable values for the bare Ising model parameters. Among the 
thermodynamical quantities explicitly calculated within this model are the internal, structural, and 
mechanical features of the DNA molecule, such as bubble correlation length and two distinct chain 
persistence lengths. The predicted variation of the mean-square-radius as a function of temperature 
leads to a coherent novel explanation for the experimentally observed thermal viscosity transition. 
Finally, the influence of the DNA strand length is studied in detail, underlining the importance 
of flnite size effects, even for DNA made of several thousand base pairs. Simple limiting formula, 
useful for analyzing experiments, are given for the fraction of broken base pairs, Ising and chain 
correlation functions, effective persistence lengths, and chain mean-square-radius, all as a function 
of temperature and DNA length. 

PACS numbers: 87.10.+e General theory and mathematical aspects, 87.15.Ya Fluctuations, 82.39.Pj Nucleic 
acids, DNA and RNA bases 

I. INTRODUCTION 

The stability of double-stranded DNA (dsDNA) at physiological temperature is due to the self-assembly of its base 
pairs: self-assembly within a same strand via base-stacking interactions between neighboring bases; and self-assembly 
of both strands via hydrogen bonds between pairs of coniplementary bases. These interactions, however, are on the 
order of magnitude of a few ksT (thermal energy) [J [3, Q and thermal fluctuations can lead, even at physiological 
temperature, to local and transitory unzipping of the double strand (see e.g. [il or the reviews @,@|). The cooperative 
opening of a sequence of consecutive base pairs leads to denaturation bubbles which are likely to play a role from a 
biological perspective, since they may participate in mechanisms such as replication, transcription or protein binding. 
For example, it has been proposed [7] that transcription start and regulatory sites could be related to DNA regions 
which have a higher probability of promoting bubbles. Indeed, the energy needed to break an adenine-thymine (A-T) 
base pair ^ AksT, connected by two hydrogen bonds, is smaller than the energy needed to break a guanine-cytosine 
(G-C) one ~ GksT (3 hydrogen bonds) [11,12 • -^^ same temperature, A-T rich sequences present a priori more 
bubbles than G-C rich ones, even though sequence effects on the occurrence of bubbles are more complex than a simple 
examination of local A-T base abundance [U, 0|- In addition to the sequence, the fraction of denaturation bubbles 
in vitro naturally depends on temperature, as well as on the ionic strength of the solution 0, In particular, the 
melting temperature T^, above which bubbles proliferate and the two strands completely separate, depends on both 
sequence and ionic strength. Another parameter that affects the melting temperature is the length of the double 
strand This is not a purely academic debate because short DNA strands (a few tens of base pairs) are involved 

in DNA chip experiments where the hybridization process is precisely affected by temperature in a way depending on 
strand length and sequence (see [ll| and references therein). 

Although the intracellular unwinding of DNA is due to active and enzymatic processes by imposing unwinding 
torsional stresses the thermally induced denaturation of purified DNA in solution has led to an intensive study 
of DNA thermal denaturation 0, [1 H, [H [11 . Mesoscopic models have been proposed to account for the thermody- 
namical properties of denaturation bubbles in DNA. The first models were Ising-like two-state models, where the base 
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pairs can be open or closed (see [J, |5| and references therein) . In the simple base-pair model, the Ising parameters 
are the base-pair chemical potential and the so-called cooperativity parameter, which accounts for the energetic cost 
of a domain wall. This type of one-dimensional Ising model is exactly soluble for homopolymers and for random 
sequences Q. More sophisticated effective Ising models, including the 10 kinds of base-pair doublets, lead to better 
agreement with experimental data [5] . Poland and Scheraga [14,] , following the work of Zimm ^ISi] , included "loop en- 
tropy" in these models, i.e. the entropic cost of closing large loops formed by destacked single-stranded DNA (ssDNA) 
in bubbles [3, . This Poland-Scheraga model is at the core of the DNA melting simulator MELTSIM algorithm 
first developed by Blake et al. [l^, [l^ ■ The intra- loop and inter- loop self-avoidance corrects the loop statistics and 
refines the sharpness of the denaturation transition [Tsl. [l9j. To get a denaturation transition in these models, an ef- 
fective temperature dependent Ising chemical potential must be inserted by hand. Recently, Peyrard et al. developed 
non-linear phonon models where the shape of the interaction potential between base pairs is more precisely taken into 
account (see [l^, [2l|,[2^ and the review [23|). By inserting estimated microscopic parameters, they have shown that 
their model can both lead to a denaturation transition, analogous to interface unbinding, and be useful for studying 
bubble dynamics. The predicted transition tem pera ture is, however, very sensitive to the model parameter values, 
and if physically reasonable values are used [2ll |24L , Tm appears to be much too high and the transition width 
much too large. 

More generally, the theoretical study of DNA denaturation can in principle be tackled on a least four levels of 
investigation determined by the amount of detail included, going from (1) quantum ah initio approaches and (2) 
classical all atom molecular dynamical simulations (26j . through (3) effective mesoscopic approaches couplin g ch ain 
conformational and base-pair degrees of freedom [13], to (4) effective statistical models for base-pairs alone (P.H,l23j 
and references therein). In theory, it is possible to move up one step in this hierarchy by integrating out the subset 
of the degrees of freedom that do not appear at the higher level, giving rise to en effective free energy at each level 
of description. Our purpose here is to show, via a minimal model for DNA homopolymers, that if one starts at the 
third level and integrates out the chain degrees of freedom, one arrives at a physically coherent level (4) explanation 
for the DNA melting transition: a bending free energy driven denaturation transition emerges naturally due to the 
entropic lowering of the energetic barrier for bubble nucleation. Working at level (3) also has the added advantage 
of allowing access to the statistics of the chain degrees of freedom (effective persistence length, mean-square-radius, 
etc.), something that obviously is not possible at level (4). 

To this end, we have recently proposed a model [2^ . which considers not only the internal coordinates in terms 
of Ising spin variables describing the open or closed states, but also external coordinates, the chain tangent vectors, 
which determine the chain configuration and depend sensitively on chain stiffness. Indeed, ssDNA is two orders of 
magnitude more flexible than dsDNA at normal salt concentration. We have shown that this difference in bending 
rigidity provides a novel explanation for the bubble mechanism formation. Further evidence for the importance 
of this bending heterogeneity include phenomena such as cyclization, loop formation, and packaging of DNA into 
nucleosomes [29l | (where denaturation bubbles facilitate bending of the otherwise rigid polymer DNA in structures 
where it coils up with curvature radii down to 10 nm, despite a persistence length of dsDNA equal to 50 nm). Our 
model incorporates precisely this dependence of the polymer bending rigidity on the state of neighboring base pairs. 
It is the discretized version of a continuous model 28] , and couples explicitly an Ising model, describing the internal 
degrees of freedom (open or closed), and a Heisenberg or discrete worm like chain (DWLC) one, accounting for 
the rotational degrees of freedom between successive monomers of the DNA chain. Its originality lies in the fact 
that the internal 2-state and external bending degrees of freedom are treated on an equal footing and therefore the 
renormalized Ising parameters obtained by integrating out the chain can be exactly calculated within the scope of the 
model. The melting temperature Tm naturally emerges and, together with the transition width AT^, are explicitly 
written as a function of the bending rigidities and strand length, which are experimentally known, and bare Ising 
parameters. In addition, the effective Ising properties (fraction of broken bases and correlation length) as well as 
the chain ones (persistence length and mean square radius) can be computed, allowing in principle direct comparison 
with experiments. An important feature of this effective Ising model is that the end monomers see an effective 
chemical potential that is lower than the interior one; this end-interior asymmetry leads in a natural way to a chain 
length dependence for T„n the transition width, and other statistical quantities. A similar model had been previously 
proposed by one of us in 2D [33], but its application to the study of denaturation bubbles in dsDNA was not made 
explicit. Recently, a similar approach has been considered in the context of the dsDNA stretching transition [3l|. 
The addition in the energy functional of the term corresponding to the external force prevents an exact solution of 
the model in Ref. [3l| and an approximate variational scheme had to be implemented. 

Beyond dsDNA or dsRNA, our coupled model can be used to describe the properties of any two-state biopolymer, 
as soon as the local bending rigidity depends on the local states. As already mentioned, the transition from B- 
to S-form of dsDNA in force experiments has been investigated in this framework [3l[. The helix-coil transition in 
poly-peptides can also be described by such a theory because the a-helix configuration is much more rigid than the 
random one [1, [l3] . 
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The present paper is a detailed account of the results summarized in a Letter [2^. In section [TTl we present the 
coupled classical Ising-Heisenberg model, which, as we show here, can be used to describe DNA thermal denaturation 
and write the partition function in terms of a transfer matrix, usual in one-dimensional statistical systems. We note 
in passing that the coupled Ising-Heisenberg model presented here displays a rich array of behavior and therefore may 
be of interest in other contexts, such ID classical spin chains or OD quantum rotators describing a diatomic molecule 
with internal states. 

Because the full transfer matrix method for the coupled model leads to relatively complex calculations, we first 
show that the model can be reduced to two effective Ising models, a path that provides a great deal of physical 
insight: indeed, these effective models allow the calculation of the free energy, as well as the Ising and chain end- 
to-end tangent-tangent correlation functions in terms of an effective Hamiltonian with temperature-dependent Ising 
parameters. A detailed solution of the two effective models is provided in section IIIII Section IIVI is devoted to the 
calculation of the Ising and chain correlation quantities, using the effective Ising models, as well as a novel correlator 
mixing both Ising and chain variables. The next sections, |V] and IVIl present in detail the full transfer matrix approach, 
which leads to the complete calculation of chain correlations and end-to-end distance; a visual interpretation of the 
expression for the 2-point correlation functions leads naturally to an analogy with a quantum diatomic molecule. Our 
theory is compared to experimental denaturation profiles of synthetic DNA in section IVlII and finite-size effects, which 
are experimentally relevant, are thoroughly examined in section IVIIIl Finally, our concluding remarks are given in 
section HXl where we also summarize our principal theoretical results of greatest interest for interpreting experiments. 
The principal symbols used in this work are defined and catalogued in Table U at the end of the paper. 



We model dsDNA as a discrete chain of N monomers (links), each monomer can be in one of two different states, 
U and B, which denote, respectively, unbroken and broken bonds. The local chain rigidity depends on the nearby 
link types. A denaturation bubble is thus formed by a consecutive sequence of B type monomers. The chain's 
conformational properties are determined by the set of N unit link tangent vectors {t^ : j = 1, . . . , N} with ||ti|| = 1. 
For simplicity, the monomer length, a, is taken to be the same for both U and B (for modeling DNA stretching 
transitions it is necessary to introduce different monomer lengths [3l|). The position of the end of the i*^ link in 
the 3D embedding space is = Xq -I- ^ji where Xq is an arbitrary starting point. The end-to-end vector 

is R = aJ^jLi ^j- The link states are denoted by the value of an Ising variable ai — ±1 (U or B) associated with 
each link. These Ising variables allow us to model a system of thermally activated defects such as the broken bonds 
that proliferate on certain macromolecules like DNA when the temperature is raised. In this case the temperature- 
dependent concentration of broken bonds is controlled by an appropriately defined chemical potential. Because we are 
interested here in the new phenomena that arise due to the coupling between the internal (Ising) and external (chain 
conformational) degrees of freedom, we will not attempt to take into account at the same time the self-avoidance of 
the chain. 

After presenting the model we explain how to calculate the partition and correlation functions connected with both 
the internal and external degrees of freedom for the coupled model. Once we have these quantities for the coupled 
Ising-chain system, we will be able to compare the results for the coupled system with those for the uncoupled one. 
We will see that the coupling can strongly modify the results, namely the average properties of the system at a given 
temperature in terms of the average concentration of closed and open bonds, average mean square chain radius and 
2-point correlation functions. 

Up to an absolute location in space a state of the chain is given by the 2A^ variables {(Tj , t^}. For a chain in 3D each 
link vector can be expressed in spherical coordinates as = (sin(6'i) cos{(j)i), sm{9i) sin((/)i), cos(0i)) and can therefore 
be defined by the azimuthal and polar angles 4>i and 9i, denoted together by the solid angle Qi = {0i, (f>i). The energy 
H[ai, ti] of a state is taken to be 



II. DISCRETE CHAIN MODEL 




(1) 



i—1 i—1 L J i—1 



The angle between two tangent vectors is given by 



cosjij — ti ■ tj = sin{6i) sm{9j) cos(0i — 0j) + cos{9i) cos{9j). 



(2) 



The first term in H is the bending energy of a DWLC with a local rigidity having the dimension of energy, that 

depends on the neighboring values of the Ising variables. We have taken, without any loss of generality, the minimum 
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FIG. 1: Illustration of the different Ising parameters appearing in the Hamiltonian (iH . Open (closed) base pairs are coded by 
a spin a — +1 ( — 1). The energies indicate the cost of opening base-pairs with respect to the ground state where all a are set 
to +1. The first line shows the cost, 2J, of a domain wall. The second line indicates the energy, 2/i, required to open a base 
pair. The third line gives the difference in stacking energy between a segment of dsDNA and a denaturated one: dark (light) 
blue cigars indicate stacked state in dsDNA (in ssDNA) and the absence of dots indicates the destacking of adjacent base pairs, 
which is already taken into account by the J contribution. 



of the bending energy to be zero independent of the values of the Ising variables. The second and third terms make 
up the energy of the Ising model [1, Hi = Hi{J, K, jl), illustrated in Fig. [TJ 

The term in J accounts for the local destacking energy (2J) of a domain wall (where Ui passes from one value 
to another). The term in K accounts for the difference in stacking energy between a segment of dsDNA and of a 
denaturation bubble. The last term gives the energy (2//) required to create a link in the state Ui = — 1 (a B link 
or broken bond). We write a dimensionless Hamiltonian l3H[ai,^li\, where /3 = l/lksT), which thus contains the 
dimensionless parameters = Pki+i^i, J = /3J, K = /SK, and /i = Pjl. The local rigidity is 




for U — U n.n. 

for B - B n.n. (3) 
for U - B or B - U n.n. 



where "U-U n.n.", etc. denotes nearest neighbor link types. In terms of the Ising field ai 

Ki+i,i = -{ku + Kb - '^KuB)<7i+i<7i + -{ku - KB)i<7i+i + cTi) + -{ku + Kb + 2kub)- (4) 

We identify the B state with two identical non-interacting single DNA strands, ssDNA, each with a local rigidity 
equal to kb/2. From Eqs. (H]) and ([U, the result of the coupling between Ising and tangent variables can already 
be predicted: both the destacking and stacking parameters J and K will be modified by the two first terms in 
Eq. (j4]) which have exactly the same functional form in cr^, whereas /2 will remain unchanged. Moreover, when all 
the bending rigidities are equal, kjj = kb = kub, the two first terms in Eq. ^ disappear and the Hamiltonian ^ 
decouples into a pure Ising Hamiltonian and a pure DWLC one (isomorphic to a ID classical Heisenberg model for 
magnetism (s^): H[ai,t.i\ — -ffDWLc[ti] + Hi[ai\. In the language of magnetism the model studied here is a classical 
coupled Heisenberg-Ising spin chain. Although pure effective Ising models have been used extensively to model helix- 
coil and denaturation (melting) transitions in macromolecules, it was necessary to introduce phenomenologically an 
effective temperature-dependent chemical potential to obtain a melting transition. A key feature of the coupled model 
used here is that a melting transition will emerge naturally in the effective Ising model obtained by integrating out 
the chain conformational degrees of freedom. 

The quantities that we will use to study the behavior of the coupled system are the mean of the internal state 
variable ("magnetization" in spin language) 



1 ^ 

c=^E^'' (5) 

i=l 

the local state average, (cti), correlation functions for the Ising variables, {<7i+rO'i), and the chain tangent vectors, 
(ti+r • ti), and the mean square radius R = (R^)^/^. These correlation functions measure the extent of cooperativity 
exhibited by the coupled system: e.g., the size of the B (U) domains below (above) the melting temperature, and the 
length scale on which the chain remains rigid. The concentration of U and B links is given by 

ipBiN^T) = =l-ipu{N,T). (6) 
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Once the type of homopolymeric DNA is chosen, the bare Ising parameters and the chain bending rigidities can be 
considered fixed, and therefore (c), Lpu, and lpb become functions of the experimental control parameters, namely 
temperature, T, and chain length, N . For a pure U chain ipu — 1 and (/3b = 0; for a pure B chain ipu = and fB = 1- 
A chain with a finite concentration of bubbles (B links) will have </?b > and the melting temperature Tm < cxd, if it 
exists, will be defined by ipu{Tm) — f>B{Tm) = 1/2 . The equilibrium statistical average of a quantity O — 0[ai,Q.i] 
that depends on the fluctuating degrees of freedom, [a-i, Vli], is given by 



p)^^ E n 

{cr,=±l} i=l 



,-pHla,.ni] 



(7) 



where 



(8) 



{<T.} ^=l 



is the partition function. The partition and correlation functions for the coupled system can be calculated using 
transfer matrix techniques. For example, the partition function can be written as 



47r 



{V\a,)(a,\P{n,,n2)W2) ■ ■ ■ {aN-i\P{nN-i,nN)\ar,){ar,\V), 



(9) 



where the transfer kernel that appears — 1 times in Eq. 



is given by 



p(a,a,+i) 



gK;[/[cos(7i-|-i_i) — IJ + J+K+p 
K[/i3[cos(7i+i,i)-l]-J 



oK!7B[cos(7i+i,i)-l]-,/ 

g[cos(7i-|-i,i)-l] + ,7-/f-/j 



It is written in the canonical base \U) = | + 1) and \B) = | — 1) of the U and B states. The end vector 

\V) ^ e^/^\U) + e-f'^^lB) 



(10) 



(11) 



enters in order to take care of the free chain boundary conditions. Different boundary conditions could be easily 
handled in a similar way, for instance for closed (open) ends, \V) = \U) (|B)) and all the following results not 
explicitly using Eq. (jlip remain valid. 

Before presenting the full transfer matrix method, we first show that the partition function and averages of any 
quantities depending only on the Ising variables can be obtained by examining the effective Ising model obtained 
by integrating over the chain conformational degrees of freedom (the link tangent vectors) in Eq. The problem 

reduces to that of an effective Ising model with an "effective free energy" H^^Jg containing renormalized parameters. 
This method works because, for the coupled Ising-chain system, the rotational symmetry is not broken (absence of 
a force term c>c • z in the Hamiltonian [1, [3l|). Hence the matrix obtained by integrating the kernel P(Jli,r2i+i) 
in Eq. © is the same for any site i. We thus are able to carry out the solid angle integrations in sequential fashion 
by using the (i + 1)*** tangent vector as the polar axis for the i*"^ solid angle integration. The solid angle integrated 
transfer matrix is 



P. 



(0) 



I. off' 



^-Ga{Ku) + J+K+n ^-~Go{kub)-J 



(12) 



where Gq{k) is the (dimensionless) Helmholtz free energy of a single joint (two-link) subsystem with rigidity n (either 
U-U, B-B, U-B or B-U): 



Go(k) 



j- = K — In 


sinh(K) 




K 



(13) 



an increasing function of k, varying linearly with k for k <^ 1 (high T) and as ln(2K) for k ^ 1 (low T entropy 
dominated regime where the spin wave approximation for the chain degrees of freedom is valid). The effective transfer 
matrix can be written in Ising form using the renormalized Ising parameters, Jo, Kq and the prefactor exp(— Fq), 



6 



all depending on chain rigidities: 



Jo = J - jiGoinu) + Go{hb) - 2GoinuB)] (15) 
Ko = K-^[Go{hu)-Go{hb)] (16) 
To = i[Go(^c/) + Go(KB) + 2Go(K[/i3)]. (17) 

In the limit of high temperature, the chain tangents will be completely decorrelated and the uninteresting renormal- 
izations of J and K arise solely from Ki+i.i (Eq. [4]), in agreement with Eqs. (fT5|l and (fT6|) . as revealed by the linear 
dependence of Go{k) on k in this limit. It is rather in the opposite limit of low temperatures and therefore strongly 
correlated chain tangents that the bending entropy driven denaturation transition arises. 
The full partition function, Z = Z^^^^, in Ising transfer matrix notation. 



4eff = (4^)"^ E (^l'^i)(^il^/,clrl^2) • • • {aN-,\P^%\aN){a^\V), (18) 



can be rewritten explicitly in terms of an effective Ising free energy, i?/ 



Z[% = (4^)^e-(^-i)i^" V e-K"i«[--l (19) 



where = Hi{Jo, Kq, i^)- 



N-l N~l N 



2 

i—1 i—1 1=1 

N-1 r- N-1 



with 



i=l i=l 

Lo = /i + Xo = M + ^-AG^^/2, (21) 

and AGq^ = Gq{ku) — Gq{kb)- Because uf^^ depends on the temperature, it cannot be considered as an Ising 
state energy, but rather as an effective free energy obtained by integrating out the chain subsystem (cf. discussion 
concerning levels of theoretical study in the Introduction) . 

When two links open the renormalized stacking energy of the links is Kq^ which is smaller than K by the difference 
in bending free energy between U-B and B-B joints. The second form of /SiJj in Eq. l(20)) shows that the effective 
chemical potential of one interior base pair is renormalized to Lq. If the gain in the one link bending free energy in 
going from U to B, AGp^'^, is greater than the intrinsic energy, 2{jl + K), needed to break a closed interior bond, 
then the effective interior joint chemical potential Lq can become negative, signaling a change in stability of U and 
B states. The end links [i = 1, N), however, feel a different chemical potential, ^ + which is larger than Lq in 

the case of interest {ku > kb)- This end- interior asymmetry, along with the extra bubble initiation energy due to the 
second domain wall, are reflected in the difference between the "effective free energy" needed to create an n-bubble 
at a chain end, 

PAGi'^l = 2 Jo - Ko + 2nio (22) 

and in the chain interior 

PAgI^^ = 4Jo + 2nLo, (23) 

which is higher than for an end link by 2Jq + Kq. As will be confirmed in section rVIIII this difference in effective free 
energy suggests that at sufficiently low T bond melting will begin at the chain ends. Written out in greater detail, 
Eq. dlSD leads to 

AG|"t^ = 4Jo + 2nkBTLo = 4 J + 2n(/i + K) - keTlGoi^iu) + Go(AtB) - 2Go{nuB)] - nkBTAG^^"". (24) 
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In the uncoupled limit {ku = kb = kub) only the first two (temperature independent) terms survive. Although 
this end-interior asymmetry plays no role in the limit of an infinite chain [N oo), it will play an essential role 
in determining how bond melting varies with temperature and bond location (melting maps) and how the melting 
temperature varies with chain size. These finite size effects are discussed in detail in section IVIIII 

Because the renormalized Ising parameters depend on the chain parameters, the coupled model does not in general 
have the same behavior as the uncoupled one. Indeed, since Go(k) is an increasing function of k and ku kb for 
dsDNA, if the difference between kjj and is sufficiently large, then at a certain temperature, T^, the effective 
interior bubble chemical potential, Lq, can vanish. Provided that this temperature be sufficiently low for thermal 
disorder to be weak and end effects due to the finite size of the chain to be unimportant, the state of the system 
will flip from nearly pure U for T < T^, where Lq > 0, to nearly pure B for T > T^^, where Lq < . Precisely 
at T = T^, there will be, on average, as many closed as open bonds and ifiB, oo — Vu.oo — 1/2, where ipB,oo{T) = 
liiiiN^oc fB{N,T), etc. This transition (strictly speaking a crossover), which can be extremely sharp under some 
circumstances (see below), can be interpreted as a melting or denaturation transition. From the above analysis, we 
see clearly that it is the difference in free energy between U-U and B-B joints, AGq^ , that drives the melting transition. 
We will see below, moreover, that if ku and kb are much greater than one, then the spin wave approximation is valid, 
and the entropy term in AGq^ dominates. We will see in section IVlIl that this is actually the case for real DNA. 

By calculating the average of the product of the cosines of the iV — 1 angles, Y[i'^os{'ji+i^i) and using the same 
technique used above for we can define another effective model, but now with partition function 

(25) 




where 



i=l I 

This expression can be written in Ising transfer matrix form as 

4off = (4^)"^ E (^l^i>(^ilA^,elrk2> • • • {ar,-i\P^^\ar,){aMV), (26) 

pSs - J -^cos{j.+,,,)P{n.,n,+,) = e-^^ ^4~K,+.h ) (27) 

Ji = J ^^[Gi{ku)+Gi{kb)-'2Gi{kub)] (28) 

Ki = i-i-]^[Gi{Ku)-Gi{KB)] (29) 

Fi ^ 1. (30) 

The function Gi{k) = — Inj/ ^ cos(0) exp[K(cos(6') — 1)]} is related to the tangent-tangent correlation function 
between two adjacent monomers (isolated 2-link sub-system) with rigidity k: 

(ti • t2)2_ii„k = (cos(6'))2_ii„k = exp [-Gi(k) + Go(k)] = u{k) = coth(K) - (31) 

which is the Langevin function [33]. It increases with k, varying as k/3 for k <^ 1 and as \~ \/k for k 3> 1. This 
asymptotic behavior corresponds to Gi(k) — Go(k) varying as ln(3/K) for k ^ 1 and as 1/k for large k. For a pure 
chain (U or B) (cos(6'))2-iink is equal to the nearest-neighbor tangent correlation function (t^+i • t^) — exp(— 1/^p) 
with persistence length 

ip{K) = -1/ \n[u{k)] = [G^{k) - Go{k)]-\ (32) 

In the high k (spin wave) approximation, we obtain ~ k, as expected. 

The partition function can be rewritten explicitly in terms of an effective Ising free energy, H^^^: 

= (47r)^e-(^-i)ri ^ g-K'-I'^'l, (33) 



where 



JV-l 



JlCTi+lCTi -I- — (Ci+i 4- Gi) 



^(al-Ka^^). (34) 
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By repeatedly using the vector identity (a • c)(b • d) = (a • d)(b • c) — (a x b) • (c x d) for a • b = t^+i • t^, etc., along 
with the property that averages of cross products, (t^+i x t^), are zero (and that |ti| = 1), the partition function 
-^i eff '^^^ written as the product of the end-end tangent-tangent correlation function, (ti • tjv), and the effective 
Ising partition function, -Zj g^: 



-"I.cff — 



AO) 

-■I,cS- 



When Ku = kb = hjjb = we recover the pure chain tangent-tangent correlation function: 



^i.cff J' exp[-G'i(K)] 
^ \exp[-Go(^)] 



(ti -t 



N 



N-l 



exp [-(Ar_i)/^p(^)] 



(35) 



(36) 



Coming back to the difference in (dimensionless) free energy AG^'^, we can show that at room temperature it is 
dominated by its entropic part. Indeed, Go{k) can be split into an average (dimensionless) energy and average entropy 
contribution, Go{k) — Eq{k) — So{K)/kB, where Eq{k) ~ f3dGo/dP = k9Go(k)/9k. Hence, the average energy of a 
two-link system can be written in terms of Gq and Gi: 



Eain) = K [1 - exp (-Gi(k) + Go(k))] = k [1 - exp {-1/^p{k))] = k [1 - u{k)] , 



(37) 



and therefore AEq^ = Eq{kij) — Eq{kb) — {ku ~ kb) — [kuu{kij) — kb'u{kb)]- Because the function u{k) tends 
to 1, AEq^ — > 0, and therefore for temperatures low enough for the spin wave approximation to be valid for both 
U-U and B-B links, we see that AGq^ ~ — AS'o^^/fcs- Indeed in this approximation, the Hamiltonian is Gaussian 
and equipartition of energy occurs: {E) ~ (3~^ ^ (E) — (3{E) ~ 1. In this case the melting transition, if it exists, 
will be driven overwhelmingly by the difference in entropy between U-U and B-B joints. As done for AG[j^^, the 
difference in (dimensionless) free energy AGf^ can be split into an average (dimensionless) energy and average entropy 
contribution using the same formulae as above. If Li = at a certain temperature, Tf°, then we can expect another 
type of "melting" transition, now driven by the free energy difference AG^^ = Gi{k,u) — Gi{kb)- We return to this 
point below. 

It should be noticed that the above processus could in principle be carried on (with increasing difficulty) to calculate 
higher order correlation functions quantities. Hence these effective Ising models give information on multi-point 
tangent-tangent correlation functions. In the next section we obtain the solutions to the two effective Ising models. 



III. SOLUTION OF THE TWO ISING MODELS 



The effective Ising partition and correlation functions can be obtained using well-known Ising transfer matrix 
techniques Q. In order to treat in parallel 2|*^|f and Z^^^f^, we introduce the index / = 0, 1 and compute the associated 
partition function -Zj^ff- This index will be useful in Section fVl where we introduce the transfer matrix approach. We 
need the eigenvectors and eigenvalues of the transfer matrices: Pj^^gl'ip'''''^) = Ai|'0^'-'), where the expressions for Pj^s 
are given in Eqs. p2l27p . The eigenvalues are 

Ai,± = e^'-^' {cosh(LO ± [smh\Li) + e"*^'] '^'} (38) 

and they obey the inequality A;.+ > A(^_. The two orthonormal eigenvectors, I'lp''-''"^) , are 

\l,+) = ^^{ai\U)+a;'\B)) and \l,^) ^ {a-'\U) - ai\B)) , (39) 

where 

7, = [sinh^iLi) + e--^''']^^^ and a; = e'^' [sinh(L,) + ji]^^^ . (40) 
The transfer matrices can be expanded in terms of the eigenvectors 

^'£ff = E vi^,^>(^,ri- (41) 
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Using the decomposition of the Ising (2 x 2) identity matrix, // = I'' ''') '''I orthonormahty of the eigenvectors 
\1,t) for each value of I [i.e., {1,t\1,t') = St-.t'], and the forms Eqs. (|18l26p for the effective Ising partition functions, 



we then find 



4off = (4^^)^ (^1 P,% \V) = {4nr E Af-i(F|?,r)2. (42) 

T = ± 



The full partition is given by Z = -Z^loff matrix elements, {V\1,t), entering Eq. can be obtained explicitly 
from Eqs. (HB), and (|4n| . 

When Kjj = Kb — kjjb, the system decouples, but the pure Ising model with temperature-independent parameters, 
J, K and /i, exhibits neither a second order phase transition at a finite temperature in zero "field" {fl + K = 0), nor 
a melting transition at a finite temperature for fl + K > 0. Indeed, there can be no melting transition because the 
inequality (pB < 1/2 holds over the whole temperature range. At low temperatures, sinh(/x + if) ^ e~^'^ , and the 
system is ordered with ipu ~ 1 and ipB — 0. As the temperature is raised, denaturation bubbles are thermally excited, 
with a cross-over when sinh(/i -|- K) ~ e~^'^, or roughly (3' = {kBT')~^ ~ {p, + K + 2J)^^. At higher temperatures 
(sinh(/i+ A') <C e^^"'), the average "magnetization" monotonously approaches a completely thermally disordered state 
with (c) — and ipu ~ tfB — 1/2. Note that this regime is not reached for DNA since as we will see below, T' ~ 8 Tm 
and the model is certainly no longer valid for such high temperatures. By contrast, the coupled Ising-chain model will 
exhibit a very different behavior, with a finite temperature melting transition. In the following we implicitly assume 
that all temperatures of interest obey T ^ T'. 

Although the eigenvectors \1,t) are orthogonal in t for the same value of /, this is not necessarily the case for different 
values of I (as we will see below, this is a consequence of the difference in rotational symmetry). In general, depending 
on the values of the temperature-dependent effective Ising parameters, Li and Ji, the eigenvectors are complicated 
mixtures of the canonical basis states, \U) and \B). These eigenvectors and their corresponding eigenvalues can, 
however, be simplified in two important limits: 

• For sufficiently low or high temperatures, below or above the transition temperature, T;°° (at which Li vanishes), 
the inequality sinh^(L;) ^ e^'*'^' is obeyed (the experimental melting temperature for infinite chains is thus 

= T^). As a consequence, the off-diagonal (domain- wall or "tunneling") terms in -P/^^g can be neglected 
and the eigenvectors reduce asymptotically to the canonical ones, with the mapping depending on the sign of 
Li: \l, +) ~ \U) and 1^, -) ~ -\B) for Li > and 1^, -) ~ \U) and 1^, +) ~ \B) for L/ < 0. In this limit of strong 
cooperativity, the eigenvalues reduce to 

Az,± ~exp(J,±|Lj|-rj) for sinh2(L,) » e"^^' (43) 
and therefore the pure U state is strongly favored if L; > and the B state if < 0, because 

^~exp(2|L,|) for sinh'(LO > e"^-^' . (44) 

By introducing the following eigenvalues for pure U and pure B 

\u = cy.^[J + ^i- Gi{ku)] (45) 
\uB = exp[J- ^ - G;(kb)] (46) 

and using the definitions of Ji, Li, and F; , these limiting forms for the eigenvalues can be further simplified 
(sinh2(Lz) »e-4'^'): 



^+-^^ for Li>0 and ^'^+"^'.5 for Li > 0. (47) 



For temperatures at or very near the transition temperature, however, the opposite inequality sinh^(i;) <C e^*^"^' 
holds and Li can be set to zero in Eqs. ([55 | - (pn|) : the eigenvectors then reduce to symmetric and anti-symmetric, 
superpositions of the canonical basis vectors, \U) and \B): 

\1,±)~^{\U)±\B)) for sinh2(L,) < e-*-^' (48) 
v2 

with eigenvalues 

Ai,± ~ e-^'ie-^' ± e"-^') for sinli^ {Li) < e"'*"''. (49) 
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FIG. 2: "Energies" ei_± — — ln(Ai,±) for Z = 0, 1 in the quantum formalism (related to the Landau-Zener problem) vs. 
temperature (for parameter values used in section IVlIl to fit experimental melting data: jl = 4.46 kj/mol, J = 9.13 kj/mol and 
K — 0). We observe that far from the two transition temperatures, the eigenvalues reduce to the limiting forms Eq. (|47|) . The 
inset is a zoom close to showing the level repulsion between the branches (0, ±). A similar level repulsion occurs near 
between the branches (1, ±). 



In this limit of weak cooperativity the eigenvalue ratio is approximately A;,+/A/.- ~ coth(J/) and the behavior 
of the system is dominated by the domain walls. The average behavior for large N, which is governed by the 
ground symmetrical state, shows vanishing average for the mean of the Ising state variable (or magnetization 
in spin language) (for I = and N oo, ipu ips — 1/2, since there is no spontaneous second order phase 
transition for the ID zero field Ising model). 

The exact results for the eigenvalues and eigenvectors interpolate smoothly between the above simplified results 
obtained far from and close to Tf° (Fig. [J) . 

Using Eq. (jlip . these limiting forms can be used to obtain simple approximations for the following important matrix 
elements : 

{V\0,+)^{ V2cosh(///2), T = T^ , (50) 
e-^^/^ T>T^ 

{V\Q,-) = { y2sinh(M/2), r = T- (51) 

T>T° 



The limiting forms for (y|0, — ) reveal that this matrix element is negative for low temperature and positive at T,"^ 
and therefore must pass through zero at a temperature T* lower than T^. This special temperature will be studied 
in detail in the section concerning finite size effects (Section IVIIip . 

If a melting transition exists at a finite temperature, T^, then Lq will go from a positive value below T,"^, through 
zero at the transition, then to a negative value above T^. This temperature dependence for Lq is similar to the 
time dependence of the uncoupled energy levels in the Landau-Zener problem, a quantum 2-state dipole system in an 
electrical field varying linearly with time [34[. It is not surprising, therefore, that the 2 branches for the "adiabatic" 
states of the Landau-Zener problem are equivalent, as the time t varies from — oo to -l-oo, to the "eigenenergies" , 
ei^± = — ln(A/^±), of the states |^,±) presented above, as the temperature varies from below the transition to above 
(with the same type of limiting behavior near and far from the transition temperature, T = Tj°°, equivalent to t = in 
the quantum Landau-Zener problem, see Fig. [5]). As in the quantum mechanics of diatomic molecules, eigenenergies 
possessing the same rotational symmetry (same I) cannot cross (level repulsion), although they do reach a point of 
closest approach at T^. As an illustration we show in Fig. [5] the Landau-Zener diagrams, e/.±(T) for I = 0,1 and 
observe level repulsion near T — and T ~ Tf° for states with same / and the possibility of level crossing for states 
with different values of I. 

From the full partition function, we define the (dimensionless) Helmholtz free energy per Ising variable of the 
coupled system, F = —N~^\nZ. The average value of the Ising state variable (or "magnetization") can then be 
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obtained for finite chains using (c) = —dF/dfj,, from which (pjj and ipB can be deduced. The expression for (c) 
simphfies in the hmit N — > oo, because only the largest of the eigenvalues, Ao,+, entering in the I — effective Ising 
partition function survives: 



N- 



dl 
9/1 



sinh(Lo) 



[sinh''(Lo) + e-4J'']i/2 



(52) 



where / = liniAr^^oo^ = — hi Ao,+. Equation (j52|) can then be used to find ipu.oc and ipB,oo- 

If Lq vanishes at a temperature, T^, low enough for the e"^''" term in the denominator to be sufficiently small, 
then the system will undergo a sharp melting transition: (c)oo will jump sharply from +1 for T < (pure U state) 
to —1 for T > (pure B state). The size of e"'*'^" term in Eq. ([52|) will determine the width of the transition region, 



AT° 



dT 



0012 



■exp[-2Jo(T° 



(53) 



which is exponentially small in Jo(T!^) when Jo(T^) ^ 1. 

In a similar manner, from the effective partition function, -Zl^^ff, we can define the free energy per Ising spin, 
i^(i) = —N^^ InZj^^jj. A quantity analogous to the average Ising "magnetization", (c), for this partition function can 



then be obtained for finite chains using 



(c(i)) 



dPW _ (c(ti.tjv)) 
dji (ti • tAf) 



(54) 



from which ip^^ = (1 



c(^)))/2 and (p^j^ can be obtained. From the definition of (c^^^), we see that it is a mixed 
correlation function that describes how the average system internal state, c, is correlated with its external state (chain 
configuration) via the end-end tangent-tangent correlation function. In the limit of an infinite chain (c^^^) reduces to 
an expression analogous to Eq. ((52|) : 



sinh(Li 



[sinh2(Li) 



-4Ji 



1/2 ' 



(55) 



where Z*-^^ — limAr^oo-f'^-'. This limiting formula is similar to the one obtained for (c) in Eq. (I52p . which implies that 
if Li vanishes at a temperature, Tf", low enough for the e~"''^i term in the denominator to be sufficiently small, then 
the system can undergo a second type of "melting" transition: (c^^^) will jump sharply from +1 for T < (pure 
c(i) "U" state) to -1 for T > Tf° (pure c^^^ "B" state). For <T < , (c^^)) ~ +1 while (c> ~ -1. This implies 
that in this temperature range (c(ti • t^v)) — ~(c)(ti • tjy). This counter-intuitive result is another manifestation of 
the non-trivial coupling between internal and external degrees of freedom. The value of the e"*''^ term determines 
again the width of the transition region: 



AT^ = 2 



9(, 



dT 



2fci3[T{-]2 



exp[-2 Ji(Tf 



(56) 



IV. ISING STATE VARIABLE-ISING AND CHAIN CORRELATION FUNCTIONS 

The average value of the local Ising variable {ui) and the 2-point {ai-\-r(Ji) correlation function can be calculated by 
starting from the expression Eq. ^ for the partition function and using the property that an insertion of a term tJj 
in the sum of products is equivalent to the insertion of the Pauli matrix in the canonical basis, 

- ( J -1 ) (57) 

at the j'^ position in the product of transfer matrices defining the partition function, Eq. (|42p . This comes from 
{a\az\(y') = crSa^a' aud the equality 

'^,{^j\pZ>,+i) = E {'^A^M){^\Pi%\^,+i) = {^M. ■ A^eirk.+i)- (58) 

cr=±l 
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We then find 



(4^) 



N 



z 



(^1 
(^1 



(0) 



I. off 



(0) 



I. off 



TV- 



(0) 



I. off 



P 



(0) 



I. off 



Pt 



(0) 



I.cff 



1^) 



By the same method used for reducing the partition function, we finaUy obtain 



(47r) 



J2 {V\0, n)Xl:;^ (0, T3|a,|0, T2)AS,,, (0, T,\a^ |o, n)x^-;~^{o, t,\V) 



(59) 
(60) 

(61) 
(62) 



with Ti = ±. The matrix elements appearing in the above expressions can be found exphcitly using Eqs. (jlip . (j39p . 
and (PU|. 

The Pauh matrix CTz, which can be interpreted as a quantum mechanical dipole moment operator, is diagonal in 
the canonical basis, \U) = | + 1) and \B) = | — 1) [see Eq. ([57)) ]. but not necessarily in the basis that diagonalizes 



P}2ff Indeed, in the / — basis we have 



^(0) 



On the one hand, the transfer matrix P/^g mixes the canonical basis states, which explains the complicated repre- 
sentation of the effective Ising partition function, Eqs. (fTS|) and (|26p. state variable average, Eq. (pTjl . and correlation 
function, Eq. (j62p in this basis. On the other hand, in the basis that diagonalizes the transfer matrix, the propagation 
between measurements is simple (no mixing), but now, in general, the "dipole" measurement process, corresponding 
to CTz, mixes such states. By directly summing Eq. (|6ip over i and using the matrix elements of Eq. (|63p . we can 
calculate {c){N,T): 



(l-(c)oo)l/^ 



(l-(c)oo)l/2 -(C)c 



(63) 



(c)(7V,T) = (c)c 



2R'v 



2Py^l-(c)^ [l_e-(^-i)/«^ 



7V[l + i?^e-(^-i)/«^] [1 



3-1/61 ' 



where 



(^|o,+) 



and ^/ is the Ising correlation length 



= l/ln(Ao,+/Ao,_),, 



(64) 



(65) 



(66) 



the typical size of minority B (U) domains below (above) T-^. Although (y|0,— ) can be positive or negative (and 
even zero), (y|0, is for physical reasons strictly positive, because both \V) and |0, -f) are linear combinations of 
the canonical basis states with strictly positive coefficients of proportionality [see Eqs. (^0]) . ([SO)) , and ([5T|) ]. The 
ratio Rv (which can therefore be negative, zero, or positive) is thus always well defined. 

The above expression, Eq. (|62p . for {ai+r'^i) can be interpreted, using "path integral" imagery, as a quantum 
mechanical measurement process over an imaginary time period of TV steps of duration 5. This interpretation is based 
on the ID classical Ising representation of the partition function of a OD quantum 2 state system [3^: the transfer 
matrix becomes the quantum propagator, P/^g ^ exp{—SH/h), where H is the quantum Hamiltonian, and the 
eigenvalues. A, of the transfer matrix are related to e , the energy eigenvalues of the Hamiltonian via e ^ — In A. In 
general the two states for the quantum system are coupled by a non-nonzero tunneling amplitude, which corresponds 
to the off-diagonal (domain- wall) terms of the transfer matrix. For instance, following Eq. ([5^ . the system is prepared 

in the initial state \V) and evolves N — r — i time steps under the dynamics determined by the propagator P/^jf , 
until a measurement of the dipole moment is performed, determined by the action of a^. The state that comes out 
of the measurement then evolves r time steps until a second measurement of the dipole moment is performed. The 
state that comes out then evolves i ^ 1 further time steps. The correlation function {at+r'^i) is thus the normalized 
amplitude that the system returns to the initial state \V) at the end of this double measurement process. 
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In the limit N oo, the results Eqs. ([6T1) - ((62)) simplify because we keep only the leading order terms (largest 
eigenvalues) which sets ti = +: 



Rv{l-{c)lY^'eM'{^-^)/^i] (67) 



where we have used {(Ji)^ = (c)oo- Using the above results we obtain the limiting form for {c){N,T) when N oo 



T. ^ ) (».-3l*.|0,..) ) (0,..|*,|0, +). (6S) 



(c)(iV,T)^^ + 2 a_vlI«S, (69) 

N^oo IM 1 — e '^-f 



In the double limit N,i ~> oo, meaning that we ignore the influence of end-monomers, expressions Eqs. (|67p and 
reduce to the simpler cyclic boundary condition forms: 



W^) „-> (c)co (70) 

{a.,+r<J,) {c)l + {l-{c)l)cM^r/^j). (71) 



■oo 



Using the limiting values obtained earlier, Eq. (I49p . we find that at the melting temperature, T, 

= -l/ln[coth(Jo)]. When g-^-'oCT-) ^ ii{T^) ~ e2"'«(^~)/2 > 1. We shah see in section EH] that ii 
can be extremely large, but finite, at T^, where it reaches its maximum value. Moving away from in both 
directions, we find that ^/ ~ l/(2|Lo|) decreases as |T — increases. When ^i{T^) ^ 1, the width of the the 
transition, Eq. ((53|l . can be rewritten as AT^ ~ ks [T^]^ /[fl£,i{T^)]. Because the system is translationally invariant 
in the limit A^,i — > oo and {{(Ji+r ~ (o'i+r))(ci — (fi))) = {(^i+rO'i) — (cj)^, Eqs. ([70| - ([7T|) show that for |(c)oc| 7^ 1 the 
Ising correlation length measures the range of correlations between spatially separated deviations of the local Ising 
spin from the average value. 

Using the same Ising model transfer matrix techniques employed above for the partition function, we can calculate 
the chain end-end tangent-tangent correlation function, (ti • tjv), which is related to the effective partitions functions 

, _ 4c^ff _ E.Af:-^(y|l,r)^ _ {V\l,+fexp[~{N-^)/e^,+]+{V\l,-fexp[-iN-l)/e,.,] 
^'■"^^ Zl% E.A^-i(V^|0,r)2 (y|o,-^)2 + (y|o,-)2exp[-(iV-l)/e/] 

where the chain persistence lengths are defined by 

Ct± = l/MAo.+/Ai,±). (73) 

This result indicates clearly that in general (ti • t^) depends on three distinct characteristic lengths: ^/ and Ci ±- 

In order to better understand the physical content of Eq. (|72p (and later results), it is useful to derive simplified 
limiting forms for the matrix elements and persistence lengths appearing therein. Using the same technique employed 
to obtain the expressions for (y|0, ±) shown in Eqs. (l50|) and ([ST]) . simple limiting forms can be derived for (y|l,±) 
in the temperature range of experimental interest, T < T^, leading to (F|l, +) ~ e^/^ and {V\l, —) ~ 

Using the limiting values obtained earlier for the eigenvalues, Eq. (|49p . we find the following limiting forms for the 
two chain persistence lengths: 

- \ 1/6 + 1/^^, X~ < r < "^^"^ T^<T<T^. ^'^> 

The limiting forms for ^ in the range < T < and ^j* _ in the range T < have a simple physical 
explanation: the effective persistence lengths, Ci±: of minority domains (B, or — , below Tm and U, or -t-, above 
Tm) tend to the typical minority domain size, ^i, when these domains behave as rigid rods 3> 6 for minority U 
domains and 3> for minority B domains). 

Is is interesting to note that the various correlation lengths can be identified with differences between the 
eigenenergies appearing in the Landau-Zener diagram (Fig. [5]): £,i = l/(eo,- — £o,+)j Ci + = l/(£i,+ ~ £0,-1-), and 
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_ = l/(ei._ — £0.+)- Hence we already observe in this diagram that ^7 reaches its maximum at which is the 
point of closest approach of the branches (0,±). In the limit of large N the expression for (ti • t^v) substantially 
simplifies and depends on only one persistence length, ^ : 

The value of N for which the limiting form Eq. ([75)1 starts to be a good approximation to Eq. ([7^ depends on the 
temperature via the weights, {V\l, +)^ and the characteristic lengths, ^/ and £,1 j-- 

V. FULL TRANSFER MATRIX APPROACH 

To calculate the general chain tangent-tangent correlation function, (t^ • t^+r), for the coupled model, we need to 
introduce the more powerful (and more abstract) transfer kernel method. This method will also shed additional light 
on the origin of the effective Ising models obtained above by first integrating out the chain degrees of freedom. To 
calculate the partition and correlation functions using this method, we need to solve a spinor eigenvalue problem in 
order to find the eigenfunctions and eigenvalues of the transfer kernel: P\'i') — A|^), or more explicitly 



(76) 



where 

I*) =*+i(f7)|t/)+*_i(f})|B). (77) 

For the pure Ising model the eigenvalues and eigenvectors can be labeled by the index r = ±. For the pure 
chain model with rigidity k and no applied stretching force (like the ID classical Heisenberg model in zero field) the 
eigenfunctions, ^pim{^) = V^Yimi^) , are proportional to the spherical harmonics, Yim{^), which are indexed by the 
integer pair (I, m), with Z = 0, 1, . . . , +00 and m — —I, . . . , +1. Furthermore, the eigenvalues for the pure chain model, 
A;, are indexed only by /, because they are degenerate in m: 



A, = 



/ 1 d \ ' rsinh(K)1 /7r\i/2 



where //+i/2(k) is the modified Bessel function of the first kind. These eigenvalues take on the values Aq = 
e^'" sinh(K)/K = exp[— G'o(k)] and Ai = \ou{k) = exp[— G'i(k)] for Z = 0, 1 and are related to the Z = 0, 1 2-link 
free energies already calculated, Eqs. (fTS]) and (pT|) [ssl ]. 

The rotational symmetry of the coupled model implies that in this case the eigenspinors can still be labeled by the 
indices {I, m; r) used for the pure Ising and pure chain models: 

(r!|*i,„,.,) =^,,„(f7)|Z,T) (79) 

and the eigenvalues, A;_t by (/,t) (degenerate in m). The eigenvalues and kets, \1,t), which are independent of the 
solid angle, fi, must be determined by solving the eigenvalue Eq. ([7S)l . In general, the eigenvalues and eigenvectors 
for the coupled system are not, however, simply direct products of the corresponding eigenvalues and eigenfunctions 
of the uncoupled Ising and chain systems. By solving the eigenvalue equation Eq. (|76p . we find that the kets, \1,t), 
appearing in the full eigenspinor, and the eigenvalues A; 7-, have already been introduced and obtained for I = and 1 
[Eqs. (|38p - (|40|) ] in the calculation of the effective Ising partition functions, Eq. (|42|) . If we define Gi{k) = — In A;, then 
the same formulas, Eqs. apply for the kets |/, r) and the eigenvalues A;,t in the general case I — 0,1, ... , +00, 

T = ±. The eigenspinors are orthonormal: 

(*i.m;r|*i',™';r') = {1,t\1',t') / —^p*^,{n)^,^{n) = 5lV 5.mm' 5rr' ■ (80) 

Once we have the eigenvalues and orthonormal eigenfunctions, we can express the transfer kernel in an abstract 
operator notation as 

+00 +1 
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and then use the orthonormahty of the eigenspinors, as well as the decomposition of unity, 

^"=E / =E E E \^l,n^;r){i'Lm;r\ (82) 

to calculate the quantities of interest in a straightforward way. As a check on the method, we can, for example, 
recalculate the partition function using the following expression: 

Z = (47r)^ W ^{V\aini){Gini\P\a2n2)---{aN-inN-i\P\<yNnN){<JNnN\V), (83) 

or in kernel product form 

Z = {Aiif {V\ P^-i \V) = (4^)^ (^l*i,™;r>'Af-i. (84) 

l.ni-.T 

Because the end vector \V) contains only the rotational ground state, |^'o,o:±) (i-C-, / = 0, m = 0), its matrix 
elements with the eigenspinors of the transfer kernel simplify to 

(*;,™;r|'^^) =<5;o'^mo(0,r|y). (85) 

Inserting this expression for the matrix element into Eq. leads immediately to the result, Eq. obtained 
previously for I — Q: 



Z = = (4^)^ ^?Kr'- (86) 



r(0) 

— 

T=± 

In order to calculate the correlation function (t^ • ti+r), we could use 

1 +^ 

t»-W = - ^ V^*„(17,+,)V'i™(r!,) (87) 

m— — 1 

which can be obtained from Eq. ([2]) and the definition of the spherical harmonics. Thanks, however, to rotational 
symmetry, the average value of • t^+r simplifies to 

(tj • ti+,.) = 3 (ti,^ • ti+r,z) = (ipioi^i) ll^lo{^i+r)) (88) 

where "010(5^) = \/3cos(0). The tangent-tangent correlation function can be written in a kernel product form similar 
to the expression for the partition function, Eq. (j83p . with the difference being that we must now insert the projection 
(or dipole) operator along the z-axis, Z — cos(0), related to ?/'io(f2) in the j — i and j = i + r positions. This operator, 
which is diagonal in the canonical basis, |crf2), has the following matrix elements: 

(o-^rJili'lcr.j+in.j+i) = -^tpio{^i)S i^i+i - Qi) ^cTi+icTi (89) 

from which the following equality can be established 

In operator product form, using Eqs. (j83l87p the correlation function then becomes 

(t, • t,+,> = 3{ATrf Z-\V\P'-^ Z P"- Z P^-''-'\V). (91) 

In the basis that diagonalizes the transfer kernel, P , the operator Z, which is not diagonal, has the following matrix 
elements: 



72_^2xl/2 /(? + l)2_„^2xl/2- 



(92) 
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which, aside from the factor (Z', t'\1, r), is the well known selection rule for quantum dipole transitions, i.e., Al = ±1 
and Am = (in, for example, the Stark effect [1^). Although {1,t'\1,t) = Stt' the matrix element {1',t'\1,t) is not 
necessarily zero for I ^ I' and r 7^ r', because in this case the matrix element is between states of different rotational 
symmetry. This result indicates that the measurement of the z-axis dipole moment can also induce a change in the 
internal state, r, of the system. 

Equation ([85|) shows that in the expression for (ti • t^+r), the Z projection operator can only connect a (^ = 0, m = 0) 
rotational state with an (^ = 1, to = 0) one, or vice- versa (as in the Stark effect for the Is state of the hydrogen atom). 
To evaluate (t^ • t^+r) using Eq. ipTj) we therefore only need one matrix element, 

(*l,0;r'|^|*0,0;r) = ^(l,T'|0,r). (93) 

By inserting the decomposition of unity between each matrix factor in Eq. (|9ip and using the orthonormality of the 
eigenspinors, we obtain 

(t,.t,;+0 = ^ E (^|0,r3)A-i(0,r3|l,r2)A^,,,(l,r2|0,ri)A^-;-(0,n|y). (94) 

When i ^ 1 and r = N — 1, using again the decomposition of unity in the |/, r) space, we recover our previous result 
for (ti -tw), Eq. 

The above expressions for (t^ -tj+r), Eqs. (PT1) - ([M)l can also be interpreted, using the "path integral" representation 
of a quantum statistical partition function, as a quantum mechanical measurement process over an imaginary time 
period of N steps. This interpretation is based on the (ID classical Ising representation) ® {ID classical Heisenberg) 
representation of the partition function of a OD quantum diatomic molecule, modeled as a 2 state rigid rotator. The 
system is prepared in an initial state \V) that is in a mixture of the internal states, t — ±, but in the spherically 
symmetric rotational ground state. This initial state evolves N — r — i time steps under the dynamics determined 
by the propagator P, until a measurement of the dipole moment along the z-axis, determined by the action of Z, is 
performed. The projections of \V) onto the rotational ground states (Z = 0, m = 0), (0,±|V^), evolve in a simple way 
because these states are associated with eigenspinors of the transfer kernel (each time step gives rise to an eigenvalue 
factor, Ao,t)- This dipole measurement excites the system from the rotational ground state to the (I — l,m — 0), 
non-spherically symmetric, first rotational excited state, along with a possible transition in the internal state of 
the diatomic molecule. The state that emerges from the measurement then evolves r time steps, until a second 
measurement of the dipole moment is performed. This measurement de-excites the system from the {I = 1, m — 0) 
excited state back down to the (Z = 0, m = 0) ground state or up to the (Z = 2, m = 0) excited state [see Eq. ([M]) ] 
again with a possible change in internal state. The state that comes out of this second measurement then evolves 
i—1 further time steps. From the representation in Eq. (|94p we see that the correlation function (t^ • ti+r) is then the 
normalized amplitude that the system returns to the initial state \V) at the end of this double dipole measurement 
process [because the final state and the (Z = 2, m = 0) excited state are orthogonal, no I = 2 matrix elements appear 
in Eq. ((94|) ]. By resolving the initial and final states, both equal to , into |0, ±) components, the sum in Eq. (jM)) 
is over all possible "time" sequences involving |0, ±). 

In order to study the limiting forms of (t^ • t^+j.), we write 



X^T ( Ao,'+ ) 



N-1 



{V\Q,' 



N-l 

, (95) 



where we have introduced the joint amplitude 

C(r3,T2,ri) = (y|0,T3)(0,T3|l,T2)(l,T2|0,ri)(0,ri|y). (96) 

When all the bending rigidities are equal to k, we recover, as expected, the pure discrete wormlike chain result, 
{ti ■ ti+r) = exp[— r/^p(K)] with given in Eq. ([32|) . 

In the limit — > 00, we keep only the leading order term cx A^^^ in Z and the surviving terms in the numerator 
2-point correlation function (ri = +) to find: 

'^O.T3 \ / Al . 



T2,T3 
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where 

^ (V|0,r3)(0,r3|l,T2)(l,T2|0,+) 

Using the effective chain persistence lengths introduced previously in Eq. ([75)1 , we can express Eq. ([97)1 in a physically 
more transparent form: 

{U-U+r) „^ exp[-r/ef,,J{C'(+,T2) + C7'(-,r2)exp[-(j-l)/6]} (99) 

T2=± 

with _ < 'Ci + s-nd the Ising correlation length already introduced in Eq. ([M)) . 

In the double limit N,i oo, the dependence on the chain ends disappears again and the expression Eq. 
simplifies to 

(t,.t,+,)^ ^(l,T2|0,+)2exp(-r/^f_,J, (100) 

T2 = ± 

which reveals the importance of the two persistence lengths, Ci ±j ^ind the two "transition probabilities" , (1, ±|0, +)^, 
for going from the ground state |0, +) to the first rotational excited state, |1, ±), with or without a change in internal 
state T. In the temperature range of experimental interest, T < rf°, (1, +|0, +)2 and (1, -|0, +)'^ can, to an excellent 
approximation, be set equal to ipu, oo and (pB,oo, respectively. When this last result is used in conjunction with the 
limiting forms for ±, Eq. ([74]) . a useful approximation is obtained for Eq. (jlOOp . valid for T < Tf°: 

{ti ■ ti+r) ^ Vu.oo exp _,_) + ifiB.oo exp (~r/^f _) . (101) 

For iV, z ^ cx) and short distances, r ^ _, wc find the limiting hnear behavior in r: 

N,l—>OQ ' 

where 



i/Cff.cF ^ (1. +io, +)ve,.+ + (1, -io, +) ver,- ^ m,oo/ei.+ + ^B,oo/e^ (los) 

is an effective persistence length for the correlation function (CF) at short distances that clearly reveals the importance 
of the shortest persistence length, under these conditions. In the triple limit N,i,r ^ oo, only one term survives: 

(t, • t,+,) ^ (1, +|0, +)2 exp {-r/e^ +) . (104) 
Exactly at T^, the above expression Eq. ()100p simphfies to 

{U-U+rhs^^ i[exp(-r/e&)+exp(-r/^P)], (105) 

which for r <^ reduces to 

(t, . u+rW:^ ^ 1 ~ ^ {i/eu + i/es) ■ (106) 

Because the inverse persistence lengths enter into Eq. (|102p . the short distance limiting behavior will tend to be 
dominated by the shortest one, above , where (1,— 10,+)^ > (1,+|0,+)^. Below T^, however, there will be 
a competition between the weights (1,— 10, +)^ < (1,+|0, +)^ and the persistence lengths, 1/Ci_ > l/^i+- 

The conditions under which these limiting expressions are valid ^proximations depends critically on the weights 
appearing in the above expressions. If, for example, (1,+|Q,+) ~ ^u,oo is sufficiently small compared with 
(1,— 10, +)^ — fB.oo at a certain temperature, then the "subdominant" term in Eq. (|100p . oc (1,— 10,+)^ (possessing 
the smaller persistence length) may actually be dominant over a wide range of r values, as shown in Fig. [3l 

Indeed, due to the coupling between bending and internal states of DNA, for realistic parameter values (cf. Sec- 
tion IVIll, the respective weights (0, +|1,±)^ associated with each correlation length change abruptly at T^: below 
T^, we have (0,+|l,+) ~ {U\U) = 1 and (0,+|l,-) ~ {U\B) = 0, thus + ~ For < T < T^, we find 
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FIG. 3: (a) Variation of the transition probabilities (1,±|0, +)^ with temperature (red solid line for + and dashed blue line 
for For T < , (1,+|0,+)^ ~ ^u,oo and (1,-|0,+)^ ~ ipB,oc. One observes that for < T < , (1,-|0,+)^ = 1 
which underlines the relevance of _ in this temperature range (with parameter values — 4.46 kj/mol, J = 9.13 kj/mol and 
K = 0, see section rVII|) . (b) Tangent-tangent correlation function given by Eq. pOOjl {N,i oo) for 3 different temperatures: 
just before the transition (dotted blue line), controlled by ^ ~ [1/^^ + 1/Ci]~^ — Ctji slightly above (solid red line), 
where the correlation length _ — is dominant for r < r* ~ 20; and after the transition {T,^ < T < T^) where the 
correlation length ^j^^ ~ disappears in favor of _ ~ (dashed green line). 
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FIG. 4: Variation with temperature of the various correlation lengths appearing in the model results: the Ising correlation 
length, 5/ (solid red line); persistence lengths of the coupled system, (appearing for long chains) and ^ (dashed dotted 
line respectively blue and green); and of the pure chains, g (dashed purple lines) (in units of a). At , the Ising correlation 
length is peaked but finite, which is related to the point of closest approach of the two (0, ±) branches (zoom of Figure [2|, and 
the effective persistence length, ^^g, rapidly crosses over from to (parameter values jl — 4.46 kJ/mol, J = 9.13 kJ/mol 
and K = 0). 




(0,+|l,+) ~ {B\U) = and (0,+|l,-) ~ {B\B) = 1 which implies ~ C^. For higher temperatures, T > T^, 
the respective weights get swapped again, but now ^^_|_ ~ ■Ci -■ These considerations lead us to introduce a critical 
distance, r*, at which the two terms in Eq. (jlOOp are equal: 

where in arriving at the last approximation we have used limiting forms that are valid when T < and assumed that 
<C Cc/j C-f (cf- Fig. 11]). When r < r* then the correlation function Eq. (jlOOp is dominated by the shortest persistence 
length, Ci Slid when r > r* the correlation function is dominated by the longest one, +. For sufficiently long 
chains and temperatures close enough to T,„, the inequality N ^ r* holds and this cross over should be clearly visible 
(see Fig. El). 
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VI. MEAN-SQUARE END-TO-END DISTANCE 

To calculate the mean-square end-to-end distance of the chain, we use the two-point correlation function obtained 
above: 

2 N N-lN-i 

= ^(t, •t,)=Ar + 2^ ^(t, •t,+,). (108) 

ij — l i—1 r—1 

When all the bending rigidities are equal to k, we recover, as expected, the pure discrete wormlike chain result: 

1 + z 2z 1 — 

R^{k) = a^NWN{u{K)) where Wn{z) = — (109) 

1 — z iV (1 — z)"^ 

with u{k) defined in Eg. pip . In the limit oo, 

More specifically, three distinct regimes can be identified: 

N, ^p(k) ^ 1 <C (freely jointed Gaussian) 

N ^p{k), 1 ^ £,p{k) ^ N (effective spin wave Gaussian) (HI) 
a^iV^, l<iV<^p(K) (rigid). 

Because Cp(k) = (.piPH) is decreasing function of temperature, the pure DWLC will go from the rigid to the effective 
Gaussian to the freely jointed Gaussian regime as the temperature is raised. 

For the coupled model the double summation in Eq. (|108p can also be carried out and we find 

R\'^ J2ti,T2.T3 C'(T'3,T2,ri)S'Ar(Ao,r3, Al,r2, Ao,ri) ( aT^" ) 

? =N + 2—^ -—W-, , (112) 



where 



z—x 



with Tn{x, y) — "^^^i^^jy^ ■ In the limit N —> oo the above complicated expression for simplifies to an effective 
Gaussian form 



R\-^ 2a' Ne,^ where ^ ^ ^(1, r|0, +)^ | + " J^'^ (114) 



is an effective " long chain" persistence length. This expression can be also obtained by using the simplified N,i oo 
limiting form for (t^ -ti+r), Eq- (llOOp . in the general formula for {R/af', Eq. (|108p . It tends to a further limiting 
form when _|_ 3> 1: 

At this expression simphfies reduces to ^^jj ~ {^fj -f ^^) /2 when ^/ 3> C^}, which is actually the case (Fig. [4]). 

For T < , the longest persistence length dominates: ~ ^\ ~ ^fj. Above T^, however, we see once 
again that there may be a competition between the persistence lengths, Ci±i ^^id the "transition probabilities", 
(1, ±|0, -|-)^, appearing in Eq. (|115p . This competition occurs now for T > T^, contrary to what was found for the 
short distance behavior of the 2-point correlation function, Eq. ()102p . because the persistence lengths themselves 
appear in Eq. ()115p . and not their inverses. For < T < Tf°, depending on the weights, the smaller persistence 
length, - J rnay actually be dominant over the larger one, +! if so, £,^g — Ci - — Cs: which is actually the case 
when we consider standard parameter values for dsDNA and ssDNA (section [VIip . as shown in Fig. 01 
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VII. APPLICATION TO SYNTHETIC DNA THERMAL DENATURATION 

Melting or thermal denaturation profiles are experimentally obtained by following the UV absorbance of a DNA 
solution while slowly increasing the sample temperature. This method allows one to follow the temperature evolution of 
the fraction of base-pairs that have been disrupted, lpb {T) . A typical profile has a sigmoid shape possibly with bumps 
that could appear depending on the DNA sequence. Different Ising-type models have been proposed 0, H, [H, [l^ for 
modeling denaturation curves by focusing on the influence of the base pair sequence, but they do not attempt to take 
into account properly the fluctuations of the DNA chains themselves. Yet, chain fluctuations increase with T and 
play a crucial role in determining melting profiles. Moreover, these fiuctuations concern both stiff' helical segments 
and flexible coils with different bending rigidities. 

In this section, we compare the model developed above, whose key element is to account for internal state ffuctuations 
on an equal footing with those of the chain, with a set of experimental data. We focus on the evolution of (^b (T) for a 
synthetic homopolynucleotide polydA-polydT. Six independent parameters appear in the theory: the polymerization 
index iV, the three Ising parameters J and ^ defined in Eq. ([T]) and Fig.[Tl and bending moduli njj for dsDNA and 
Kb for ssDNA. Note that we have also introduced a bending rigidity kub for domain walls. However kub appears in 
the theory only in the effective cooperativity parameter Jg. Thus changing kub is equivalent to varying the bare J, 
i.e. the energetic penalty to create a wall. Without any lost of generality, we choose to fix kub = k,u- Moreover, we 
choose free boundary conditions for the end monomers, which is valid unless their state is fixed by the experimental 
conditions [STj (although any type of boundary conditions can be treated using our model). Of the six parameters, 
three are determined experimentally: N , ku, and kb- Moreover, there is evidence that stacking interactions in dsDNA 
and ssDNA are of the same magnitude which justifies the choice oi K = adopted below 38]. 

Figure [5] shows ipB{T) for a polydA-polydT of molecular weight = 1180 kDa in a solution of 0.1 SSC (0.015 M 
NaCl + 0.0015 M sodium citrate, pH 7.0) taken from 
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FIG. 5: Fraction of broken base-pairs for a polydA-polydT vs. temperature (solution of 0.1 SSC, A'' = 1815). The solid line 
represents the theoretical law for fi — l.GAkBTm and J = 3.35 ksTm where Tm ~ 326.4 K. The case N oo for the same 
parameter values corresponds to the broken line. 

In order to compare the data with our model predictions, we choose the experimental values persistence lengths, 
~ 50 nm and iP, ~ 1 nm at 300 K, which lead to ku = ids/a = 147 and kb = 2£P,/a = 5.54 at T = (taking 
a = 0.34 nm for one base-pair size and a factor of 2 for two flexible segments in parallel per coil segment). The two 
remaining parameters fl and J are determined by fitting the experimental data. The solid line in Fig. [5] corresponds 
to /i = IMkBTm ^ 4.46 kJ/mol and J = 3.35 kBTm ^ 9.13 kJ/mol leading to Tm = 326.4 K. We can then deduce 
several thermodynamical features. Noting that the bare enthalpy for creating one A-T link is, in our model, 2/i, we 
find a value very close to the experimental value of 10.5 kJ/mol Although the value of J is more difficult to 
interpret, our result J ~ 2/i is consistent with the idea that stacking interactions make the dominant contribution 
to DNA stability Chain fluctuations do not only renormalize the effective free energy, 2Lo, required to break an 
interior base-pair, but also the cooperativity parameter Jq: the latter varies almost linearly with T following Eq. (jl5p 
contrary to previous theories where J was taken as constant and supposed to be purely enthalpic in character [J] . We 
have for the total cooperativity parameter Jq = 4.17 kBTm at T = T^, which shows that the bending contribution is 
roughly 25% (remembering that we have chosen kub = ku). The model flt thus leads to parameter values in accord 
with experiment. Our model predictions for experimentally accessible A-T pair quantities are also in agreement with 
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accepted values 0, H^: i) the loop initiation, factor, a^i = e"^"^" ~ 10~^ at T^; and at physiological temperature, 
Tph, ii) an interior single base-pair opening probability (/^^(Tpii) ~ 10~® with a bubble initiation barrier of llksT , 
and iii) a free energy of O.lSkBT for breaking an additional base-pair in an already existing bubble. In reality, the 
fitted values of fi and J implicitly compensate for effects like loop entropy explicitly left out of the model [4]. As 
shown in Fig. 9 of effective Ising models without loop entropy, like ours, can be considered to account implicitly 
(and approximately) for loop entropy, provided that one allows for loop entropy contributions to both J and K. This 
loop entropy renormalization of the Ising model parameters will depend on the value of the loop entropy exponent k 
and the chain length N and could allow for a simple app roximate way of accounting for the influence of loop entropy 
within the framework of an effective Ising model (cf. [40|). This renormalization probably explains why our the model 
value for (Tli is at the low end of the accepted spectrum. 

In Fig. the curve corresponding to the thermodynamic limit {N — s- cx)) is shown for the same parameter values. 
In this case, (Pb{T) is given by Eq. ([5^ and the value of the melting temperature is obtained analytically as a function 
of T,^ using Lo{T^) = which is given in the limit of low temperature {kb ^ ksT) by 

fcBr-^2 (116) 

\n{Ku/KBj 

Hence, the melting temperature is reached when the enthalpy required to create a link is perfectly balanced by 
the difference in (entropy dominated) free energy between the two types of semi-flexible chains (U or B). Another 
quantity which has an experimental relevance is the width of the transition. In the thermodynamic limit this width 
is narrow, but nonzero, due to the large but finite cooperativity parameter: AT^ cx i/^i{T^) ~ 2 exp[— 2 Jo(r^)] 
[see Eq. (j53p ]. Hence, the thermodynamic limit clarifies the role of the two free model parameters: in conjunction 
with the experimentally known bending rigidities, /x sets the melting temperature and J fixes the transition width. 
This is in contrast to previous Ising- like models, where three fitting parameters were used J, dL^/dT, and with 
Lq assumed to by a linear function of T Q • 

Within the scope of our model the measured transition width is indicative of a very long Ising correlation length, 
^7, near the transition temperature, much larger than the pure U and B persistence lengths; therefore typical helix 
(U) and bubble (B) domains (of size ^ ^/) are flexible within a small temperature window near the transition. 

Included in the predictions of our theory are mechanical and structural features of the composed chain, such as 
persistence length or mea n sq uare end-to-end radius, R. This differs from purely Ising-type models [4, 5] and non- 
linear microscopic models [20l[2]| where only thermodynamical quantities related to base-pairing are available. The 
variation of the effective persistence length (and thus the radius of gyration for long chains) vs. T is shown in 
Fig. m It varies from for T < to for T > T^. Since the transition is very abrupt, we suggest that the 
denaturation transition can also be followed experimentally by measuring directly the radius of gyration, for instance 
by tethered particle motion (4l| . light scattering, or viscosity experiments. For instance, since the relative viscosity is 
proportional to cdna^^'^ (where cdna is the DNA concentration) , it should clearly exhibit an abrupt thermal transition 
at a given cdna and N. Such a transition has indeed been observed for the viscosity of synthetic homopolynucleotide 
solutions [13 , in qualitative agreement with Fig. 21 

In fitting our model to experiment for chains of length N = 1815, we have found that finite size effects play an 
important role (Fig. [5]). In the following section, we investigate such effects in detail. 



VIII. FINITE SIZE EFFECTS 

It has been shown experimentally that DNA thermal denaturation varies with chain length, N [l3| . In this section, 
we carefully study the effect of chain ends on the denaturation transition. In Fig. [5] are shown denaturation profiles 
for various chain lengths from N — 100 to — > oo and the fixed parameters values p, = 4.46 kJ/mol, J = 9.13 kJ/mol 
and K = used in the previous section to fit the melting data. Within the scope of our model one observes that 
i) the melting temperature Tm{N) is a decreasing function of N, varying as [Tm{N) — T!^]/T^ ~ 1/{N — 1); ii) all 
the denaturation curves intersect at a temperature T* at which ipB — 0.03; iii) the transition width ATm(Af) is a 
decreasing function of A^. 

Concerning points (i) and (ii), the observed behavior for the coupled system with the present parameter values is 
directly related to the model result that T* < T^, which is radically different from the behavior found for the simple 
Ising model (for which melting curves, ipB, are strictly decreasing functions of A'^, because formally, T* — oo when 
Ku = Kb = kub)- The present behavior for melting maps is also very different from the predictions of older empirical 
Ising-like models of denaturation and Helix-Coil like transitions @, for which the chemical potential /i appearing in 
the end vector \V) is incorrectly identified with Lq. This identification results in melting curves independent of TV, 
i.e., Tm = T*. 
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FIG. 6: Linear-Log plot of melting curves for A'^ = 100 (dashed dotted green line), 500 (dotted purple line), 1815 (solid blue line), 
oo (dashed red line) in decreasing order at low temperature, T < 326 K (parameter values jl = 4.46 kj/mol, J — 9.13 kj/mol 
and K — 0). We note that all the curves intersect at T*, as discussed in the text. Tm is defined by ifs ~ 0.5. Inset: Log-Log 
plot of model results for the shift in transition width ATm — AT^ vs. polymer length. Dots correspond to the model results 
and the solid line is a law in 1/A^. 




Concerning point (iii), the transition width roughly follows the law (AT^ — AT,^) ~ XjiN — 1), which is a classical 
result for finite size systems where fluctuations decrease in the thermodynamic limit. One observes that even for 
a long polymer, iV ~ 10'^, finite size effects are important. For very short chains, e.g., N < 100, such effects get 
amplified and we predict a transition width as large as 50 K for N = 10. This point is crucial, since it has been 
observed experimentally that for polydA-polydT inserts between more stable G-C rich domains, melting curves are 
much wider for very short DNA chains {N ^ 10 bp) [s^l with a width that decreases with increasing N (observed for 
60 < < 140 in [lOl)- In such experiments, the nature of end monomers clearly becomes extremely important. 

For a given N and T the local site dependent bubble opening probability, or melting map, 

^B. = (117) 

can be obtained from (ai) given in Eq. (j61[) . In Fig. [7] (/Ss^i is plotted for six different temperatures using the same 
model parameter values employed in Fig. [Sj we observe that below T* the chain unwinds from the ends, whereas 
above this temperature an interior bond has a higher probability of being open than an end one. Far enough below 
T* the melting curve heals rapidly to a plateau value close to ipB,oo on a length scale on the order of ^7 <^ N. At 
physiological temperature, 310 K, the interior bond opening probability is ~ 10~^ in agreement with the experimental 
value for long runs of A-T pairs (which is an order of magnitude lower than randomly placed A-T pairs) [3|. At this 
temperature the end bonds have opening probabilities two orders of magnitude greater than the interior ones. 

At T* the melting curve is perfectly flat. This result, which is independent of chain length N, indicates that 
each Ising variable can be considered to fluctuate independently: (ai) is constant independent of i, despite a 2- 
point correlation function that does not factorize [{aiai+r) 7^ {'^i){'^i+r), cf. Eq. (|7ip ] and a large Ising correlation 
length (^/ ^ 1, see Fig. U). Indeed, the influence of the renormalized stacking energy (~ Kq < 0), which favors 
bubble formation, exactly compensates that due to the renormalized destacking (~ Jg), which suppresses bubble 
formation, and therefore (pB,i{T*) = [1 — tanh(/i/(fcBr*)]/2 ~ 0.03, which results from taking = 1 or taking 
Kq = Jq ~ for arbitrary N. In some ways this compensation leads to an effective non-interacting Ising system with 
{d{c) / dN)T- — {dipB/dN)T* = 0, which explains why the melting curves for different values of N cross at T* in 
Fig. [6l Below T* the combined effects of the renormalized stacking energy and entropy gain favoring interior bubbles 
are too small to overcome the destacking energy cost associated with an extra domain wall and the chain ends unwind 
first. Since Kq becomes more negative with increasing T faster than Jp increases, a temperature T* is reached where 
the renormalized stacking and destacking effects just compensate. Above T* the situation is reversed and the opening 
probability is higher in the chain interior. For arbitrary N and T, the thermodynamic chemical potential, defined by 
ft = {dF / dN)T can be related to (c) via Maxwell-type relations, leading to {dfi/ djl)T.N = N{d{c) /dN)T ~ (c). At T* 
this general relation simplifies to {dp,/dfi)T* ,n = — (c)|t' = {d f / dji)T' ,n , characteristic of a non-interacting system. 

Upon examination of Eqs. (pT|) and (|86p. we see that T* is determined by the condition that {ai) ~ {c)oo, which 
is obtained when the end vector \V) is identical to the eigenket |0,-f) and orthogonal to |0, — ): (F|0,-|-) — 1 and 
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FIG. 7: Average melting maps for different temperatures for A'^ — 1815 and the same parameter values used in Fig. (5] plot 
of the fraction of broken bases, '^B,i as a function of the base position i for, in increasing order, 310 K, 0.99 T* = 322.87 K, 
T* = 326.13 K, = 326.24 K, r„ = 326.4 K and 1.01 T„ = 329.66 K. 



(V|0,— ) — 0. Physically, this means that the coupled Ising-chain system is in a pure state, |0,+), that mixes the 
canonical states in a special way. The temperature T* can be obtained by solving (V^|0, — ) — 0. Using Eqs. pT|) 
and p9l) - (|40)) this translates into e'^ = e~^'^°{smh{Lo) + [sinh^(Lo) + e^^'-'"]^/^}. After some manipulation using 
Eq. dig), this can be shown to be identical to {c)^{T*) = {c){N = 1,T*) = tanh[/i/(fci3T*)]. Furthermore, Eqs. §^ 
and ()94|) show that the Ising and chain 2-point correlation functions, (<Ti • (Ti+r) and (t^ • t^+r), also get simplified at 
T*: the approximate forms, Eqs. ([7T|) and (jlOOp valid in general only for N,i oo, become exact for arbitrary N 
and i at this special temperature. It is clear that at T* the coupled system behaves as if there are no end effects and 
finite chains have the same behavior as an infinite one. 

To shed additional light on this mechanism and illustrate the important role of internal bubble entropy for long 
chains, we now study an infinite chain and compare ipB,int — hmjv,i^oo V's.i with (ps.cnd = ^^'^n—^oo'^b.x- Equa- 
tion ((67)) shows that ^/ plays here the role of a healing length, over which end effects relax (see Fig. [7]) . The ratio of 
matrix elements appearing in Eq. (j67p gets simplified in the following way for special values of T: 



which shows that for very long chains (TV > at T^{> T*) ipB,cnd = [1 - tanh[/i/(2fci3r^)]/2 ~ 0.16 < 
'fB,oo{Tm) = 1/2, revealing an internal opening probability more than three times higher than an end one. In 
Figs, m and m however, we observe that for T = Tm and , > N — 1815 (cf. Fig. [4]), and therefore end effects do 
not get damped out near the center of this finite chain. Indeed, at the opening probability fB,i near the middle 
of a chain of length A'' = 1815 is much less than the value of 1/2 holding for an infinite one. For A^ = 1815 we still 
observe, however, noticeable differences (~ 10 to 20%) between interior and end opening probabilities near Tm- 

The one-sequence-approximation has been defined by Poland and Scheraga [H and consists in neglecting the many 
small bubbles which eventually collapse and considering only one large thermally excited bubble. It is valid for 
temperatures sufficiently far below T,"^. For N ^ oo, ips, int and ifB, end can be estimated in this approximation by 
summing over, respectively, all the interior and end bubbles containing the fixed site in question. In the interior case 
we thus have 




(118) 



At T^, Eq. ^ simplifies to 




(119) 



OO OO 



g-4Jo 



VB,int — 



^nexp -/3AG|:) = e'^''" ^ n e 



-2nL, 



— 



4sinh^(Lo) 



(120) 



n— 1 n=l 
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FIG. 8: Zoom of the melting map shown in Fig. [7] for (lower curve) and Tm (upper curve). 

The factor of n in the sum is entropic in nature and equal to the number of ways of placing a fixed interior site within 
an interior n-bubble. In the end case, where there is no entropic factor, 

^ r , M ^ „-(2Jo+Ai) 

. J: exp -^AG« = Z e--. = (121) 

71=1 n=l ^ ^' 

It is important to note that for the model parameters employed, the additional dimensionless free energy for breaking 
an additional base-pair in an already existing bubble, 2Lq ~ 0.18 at physiological temperature, is less than one and 
much smaller than the bubble initiation energy cost for an end bubble, ~ 2 Jq ~ 8 (and a fortiori for an interior one, 
~ 4 Jo). This implies that even at physiological temperature, where the probability of bond opening is very small, 
bubbles covering a wide range of sizes contribute to the sums in Eqs. (|120p - (|12ip : both (rt)int — coth(Lo) — 1/io — 6 
and (n)ond = e^^° /[2smh{Lo)] ~ l/(2|Lo|) ~ 3 are of the same order of magnitude as which varies as l/(2|Lo|) 
when |Lo| ^ 1- Because Lq decreases with T (going to at T^), both (pB.int and (pB,end are increasing functions of 
temperature. For x < T* the bubble initiation energy cost dominates and fB,int < fB,cnd- Thanks to the entropic 
factor, however, ipB;nrt increases more rapidly than ipB,cnd and the two curves cross over at T*, an estimate of which 
can be obtained by equating the above two one-sequence-approximations for both these quantities. As a check on the 
preceding discussion, the one-sequence-approximations for infinite chains obtained above for end and interior opening 
probabilities can be shown to be in exact agreement with what one gets from the definition of ipB,i, Eq. Ijll7[) . and 
Eq. ([57)1 by using the low temperature approximations for {V\0, —) / {V\0, +) ~ — e~'^ [Eq. (jllSp ] and (c)oo, Eq. ([5^ . 
i.e. expanding to lowest order in e~^'^° , assuming that e^'*''" <C sinh^(Lo), which is the formal criterion for the validity 
of the one-sequence-approximation. 

Using the one-sequence-approximation, it is now easy to see how the cost in l oop entropy associated with internal 
bubbles will modify the above results. In the so-called loop entropy models [1, [ij, [l3| the internal bubbles formed 
by the single strands are visualized as one polymer loop, whose entropic cost has been estimated first by Zimm [l5{ . 
Hence, although ipB, end does not change, the interior opening probability does, becoming 

oo 

V>¥^int ^ E '^/("o + '^)' exp[-/?AGL",)], (122) 

n=l 

where we have adopted a common simplified form for the loop entropy factor, /le('^) — {no + rt)^'° parametrized by a 
constant ng, which may be as large as 100 [l3|, and an exponent k, usually assumed to be in the range 3/2 < k < 2.1, 
depending on the extent to which chain self-avoidance is taken into account [Tsj . A large value for hq would severely 
reduce the importance of loop entropy for short chains and probably reflects the presence of strong bending rigidity 
effect (an important open question concerns how to incorporate bending rigidity into /le in a physically correct 
way). Loop entropy clearly lowers the probability of interior bubble opening and will lead to an increase in T*. 
The situation is further complicated if strand sliding, which may be important for periodic DNA, is taken into 
account. For homopolymeric DNA like polydA-polydT, strand sliding leads to a modified loop entropy exponent, 
k' = k — 1 0, [3, 13, resulting in a significant decrease in the importance of loop entropy. By the way, it also 
minimizes the importance for homopolymeric DNA of recent claims that a true first-order phase transition should 
occur for infinite chains because k appears to be greater than 2 when self-avoidance is fully taken into account p^ . 
The combined effects of loop entropy and strand sliding will lead to increases in both T* and T^. Neither the 
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theoretical nor the experimental situation concerning Tm{N) is entirely clear for homopolymeric DNA with free ends 
and further careful experiments are clearly called for. If loop entropy and strand sliding are included in the coupled 
Ising-chain model, T* might become higher than T^, which would imply that Tm{N) would increase with N 3], 
unlike what what we find to occur when these two effects are neglected. 

In the future, we intend to examine these questions by incorporating loop entropy and strand sliding directly into 
our DNA model. With all other factors being equal, adding loop entropy and strand sliding increases the stability of 
the closed state, resulting in sharper melting curves and higher values of T^{N) (see Figs. 9 and 10 of fi| and Soi). 
The importance of this loop entropy contribution would change with chain length iV and may lead to an increase in 
Tm{N) with increasing N [lj| at least for a certain range of chain sizes . 



IX. SUMMARY AND CONCLUSION 

This paper presents a novel theoretical model of DNA denaturation, already introduced in , which focuses on 
the coupling between the base-pair link state (unbroken or broken) and the rotational degrees of freedom of the 
semi-flexible chain. The Hamiltonian includes local chain bending rigidities whose values depend on neighboring 
base-pair states: around 5 fc^T for bubbles and 150 fcsT for connected base-pair segments. Because of the rotational 
symmetry, the model can be rewritten in terms of an effective Ising Hamiltonian by integrating out the rotational 
degrees of freedom of the chain. Hence, our model yields considerable insight into the empirical temperature-dependent 
parameters used in previous Ising-like models [4j|. In particular, the melting temperature r„i is no longer a fitting 
parameter, but emerges naturally as a function of: (i) experimentally known bending rigidities, kjj and kb', (ii) the 
bare energy required to open a base-pair, 2fl; (iii) the bare energy of a domain wall, or destacking, 2J; (iv) the 
difference in bare stacking energy between ss and dsDNA, 2K] and (v) the polymerization index, N. Moreover, our 
model allows structural features of the DNA chain, such as the mean size R, to be calculated as a function of T. An 
abrupt transition for R is found at Tm and explains, at least qualitatively, the thermal transition observed in viscosity 
measurements (43 |. 

From an experimental perspective, our results obtained from exactly solving the coupled model can be summarized 
as follows. First of all, we propose formulae for chain free-boundary conditions, this information being encoded in 
the end vector \V) given in Eq. (|lip . However, any other boundary condition can be treated following the same 
route, even though we shall not detail the calculations here. For example, a polydA-polydT sequence of length N 
sandwiched between more stable G-C sequences [s^ can be seen near its melting transition as a DNA of length N, 
with fixed boundary conditions, resulting in an end vector \V) = \U). 

Once boundary conditions are set, our model predicts melting profiles, as measured for example from UV absorbance 
experiments. Even if the model relies upon six microscopic parameters, as discussed in section VII, most of them are 
known experimentally, including the strand length N, and only two of them must be extracted from melting profiles: 
/2, the bare half-energy required to break a base pair (which can also be estimated experimentally and J, the 
cooperativity parameter that indicates the cost of creating a domain wall between unbroken and broken base pairs 
(recent experiments on single DNA molecules aim at determining this quantity, see Melting profiles, giving the 

fraction of broken base pairs, fs, as a function of the temperature T, are determined through the average of the Ising 
state variable, (c), because ipB{N,T) = [1 - (c)(A^, T)]/2. 

For infinite chains, {N — > oo), the expression for (c)oo(T') is rather simple (Eq. (|52p ): 

/ \ (rr\ sinh(£o) 

^ [sinh2(Lo) + e-4Jo]i/2 ' (123) 

where the renormalized parameters Jqj Kq, and Lq = /i -I- Kq are given in Eqs. (fT5|) . (fT6|) . and (|2ip . In these 
latter equations, the function Go(k) has a simple algebraic expression (Eq. (|13p ) that reduces to a purely entropic 
contribution, Go(k) — ln(2K), in the physically relevant low-T (spin wave) approximation. From the melting profile 
fB,oo{T), the melting temperature, T^, is defined by fs^oo = 1/2, in other words by Lq — 0. The finite transition 
width is estimated by Eq. AT,^ ~ 2 ^^[r^]^ exp[-2 Jo(T^)]//i. 

For finite length strands {N finite), the expression for (c)(A^, T), Eq. (I64p . is more complex, because N has an 
infiuence on ipBiN,T) and thus Tm- In addition, the interplay with mechanisms, such as loop entropy, not taken into 
account in this study is nontrivial, as discussed in detail in section IVIIII At the level tackled in the present paper, we 
obtain a simplified expression for (c) {N, T) when TV is large, Eq. ([S^ . which simplifies even further when iV ^ ^/ ^ 1 
(Fig. [71): 

{c){N,T) {c)^+2Rv{^j/N) {1 ~ {c)l)^/\ (124) 
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where Rv, Eq. (j65p . which is a ratio of matrix elements pertaining to end effects, simphfies at certain special temper- 
atures, see Eq. (|118p . Note that finite-size effects are still important for sizes of several thousands of base pairs (see 
sections IVIII and IVIIII and Figs [7] and [8]) and are not a purely academic debate. 

Three important correlation lengths can be calculated in the framework of our model. On the one hand, the Ising 
correlation length gives access to the typical size of bubbles in the low temperature regime (T < T„i), as well as 
to the typical size of unbound regions for T > Tm- This quantity is calculated in Eq. ([66]) and assumes a simplified 
form at T^: ^i{T^) ~ exp[2 Jo(T^)]/2 > 1, when Jq{T^) > 1. On the other hand, one effective chain persistence 
length, ^gjj (-,p ~ \'fiu,ool^ij + ¥'B.oo/^s]"^' provides information on the short distant behavior of the chain tangent- 
tangent correlation function, Eq. (jl02p . and the other, ^^^j, provides information on the typical chain conformations, 
in particular its mean-square-radius: 

(R2> ~ 2a'Ne,s - la^N^^uMu + ^b,ooC^), (125) 

where a is the monomer length (0.34 nm) and the approximation is valid for very long chains (large N). Kno wing this 
quantity is of primary importance when interpreting data from Tethered Particle or Tweezer experiments (4ll |45| |46| , 
Atomic Force Microscopy 47] or viscosity measurements t42i]. The results that we have obtained here for the chain 
tangent-tangent correlation function and mean-square-radius are very different from what one obtains by solving 
a quenched random rigidity model, where the local joint rigidity can take on one of two values, ki and K2, with 
probability ipi and ip2 — \ ^i- In this case there is only one effective correlation length: (t.^ • t^+r) — e~''/'>r'an and 
(R2) ~ 2a^N£P^^, where = -1/ ln[(^ie-i/«p(''i) -I- (^ze'^/^''^''^'] and ^p(k) is given by Eq. ([32l). 

The foregoing analysis makes allowance for neither solvent entropic (hydrophobic), nor electrostatic effects, which 
might not only change the actual value of the bare Ising parameters J, K and /i, but also lead to additional entropic 
contributions. To go further concerning solvent entropic effects would require molecular dynamics simulations, which 
is beyond the scope of the present approach. The electrostatic effects in DNA melting are two-fold: an entropic 
contribution arising from the difference in electrostatic energy between states U and B and an an enthalpic contribution 
arising from counterion release. At equilibrium the two effects partially compensate and the remaining contribution 
is, using a simple thermodynamic approach in the low salt limit [isl l49l|. approximately equal to AGci = kBT£h{TB — 
T[/) ln(///o), where / is the ionic strength (/q = 1 M), £b = / {AnekBT) the Bjerrum length (e is the elementary 
charge and e the water dielectric permittivity) and tu and tb are the linear charge densities of the pure U and B 
chains, respectively. Although this is an approximate result and the determination of tb — tjj remains controversial, 
these two contributions should be included in a more refined model. A natural extension emerging from this study 
concerns the ionic strength dependence of DNA melting profiles and effective persistence lengths. 

The current rapid development of force experiments in magnetic or optical tweezer traps |45i] , "Tethered Particle 
Motion" experiments [4l| or even more recently atomic force microscopy ones (47j , presents a formidable opportunity 
to investigate directly the elastic properties of DNA strands as a function of temperature, salt concentration and 
length TV. Indeed, despite the pioneering work by Blake and Delcourt [l0|, a systematic experimental study of the 
effect of N while controlling the nature of chain ends, is lacking, especially for homopolymers. This would be a way 
to discriminate between the different models and also shed light on the role of loop entropy, which is neglected in the 
present model. 

The approach developed here can be extended to bubble dynamics. Very recent work [s^, [Ell studied the 
growth of already nucleated bubbles using the Fokker-Planck equation applied to the Poland-Scheraga model (i.e., 
an effective Ising model including loop entro py) . The agreement with experimental results obtained by fluorescence 
correlation spectroscopy is remarkably good [39]. The issue of bubble nucleation is, however, not solved and will be 
continued to be explored in the near future. In addition the mutual influence of thermally excited bubbles and chain 
flexibility should play an important role in determining global DNA conformations and strongly influence looping 
dynamics [46l. [s^. 
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Symbol 


Quantity 


Mathematical definition 


Reference 


T 


melting temperature 


ipuiTm) = ^B{T^) = 1/2 


Sees. U M 


iV 


chain length 




Sec. HI] 


U/B 


unbroken/broken base pair 




Sec. HI] 


'PU,B 


Iractions of U's and B's 




Eq. © 




internal degree of freedom (Ising variable) 


±1 


SecHI] 




chain unit tangent vector (Heisenberg variable) 


llt.ll = l 


Sec.m] 




local chain bending rigidity (coupling) 


K(7, Kb or K(7_B 


SecHI 


J 


half-energy of a domain wall 




SecHI] 


A" 


difference in stacking energy between ds and ssDNA 




SecHI 




half-energy required to open a base-pair 




SecHI 


K, J,K, n, . . . 


adimensional energies (in units of fcsT) 


K — I3k, . . . 


SecHI 


GoW 


free energy of a single joint of rigidity At 


hi — ln[sinh(K) / k\ 


SecHI 


Jo 


renormalized (effective) Ising parameter J 


J ~ \[Go{ku) + Go{kb) - 2Go(kc/s)] 


Eq. ^ 


A'o 


renormalized (effective) Ising parameters K 


K - \[Go{nu) - Go{kb)] 


Eq. CSJl 




effective chemical potential to open an interior base pair 


^J. + Ko 


SecHI 




infinite size Ising "magnetization" 


sinh(Lo)/[sinh2(Lo)+e-*-^«]i/^ 


Eq. dSUl 


fU,oo, 'PB,oo 


fraction of unbroken (U) and broken (B) bonds (A'^ oo) 


[l±(c>oo]/2 


SecHU 




site i bonding opening probability (melting map) 


[1 - {ai)]/2 


Eq. (fTTfll 




infinite size (A'^ — > cxd) melting temperature 


Lo{T^) = 


SecHI 


AT°° 


transition width (A'^ — + oo) 


2\d{c)oo/dT\-^ at & Figs. [5] and |6] 


Eq. ((53]l 




Ising correlation function 




SecHVl 




Ising correlation length 


e^-'" /2 at Sz Fig. g] 


SecHVl 




chain correlation function 


Fig.Hb) 


SeclV] 




persistence lengths 


Fig.m 


Eq. (fTOOll 


Scff 


effective persistence length 


Fig.S] 


secrvn 


R 


chain mean square (or gyration) radius 


(R2)1/2 


Sees. HTTV^ 




crossover temperature for finite chains 


d<PB/dN\T- = 


Sgc.IVIIII 



TABLE I: Index of the main symbols used throughout the paper, with their mathematical definition, and reference. 



